\documentclass[11pt]{article}
%\documentclass[11pt,draft]{article} 
%\pdfoutput=1 
%\usepackage{showkeys}
\usepackage[utf8x]{inputenc}
%\usepackage[french]{babel}
\usepackage[english]{babel}
\usepackage{graphicx}           
%\usepackage{url}
%\usepackage{eurosym}
%\usepackage{makeidx}
\usepackage[plain]{algorithm}
\usepackage{algorithmic} 
%\usepackage{color} 
%\usepackage[margin=3.3cm]{geometry}
%\usepackage{setspace}
%\usepackage{datetime}
%\usepackage{pifont}
%\usepackage{tikz}
\usepackage[numbers]{natbib}
%\usepackage{endfloat}

\usepackage{amsmath,amsfonts,amssymb,latexsym}
%stmaryrd}

%\doublespacing



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%---------------------------------------------------------------------------
\input lib/e-macros
\input lib/macros
%\input lib/algorithmic
\newenvironment{psmallmatrix}{\left(\begin{smallmatrix}}{\end{smallmatrix}\right)}
\newcommand{\Sin}{S_{\textrm{\tiny in}}}
\newcommand{\mumax}{\mu_{\textrm{\tiny max}}}
\newcommand{\Mfilter}{M_{\textrm{\tiny filter}}}
\newcommand{\deltafilter}{\delta_{\textrm{\tiny filter}}}
\newcommand{\deltasimu}{\delta_{\textrm{\tiny simu}}}
\newcommand{\Msimu}{M_{\textrm{\tiny simu}}}
\newcommand{\nmax}{n_{\textrm{\tiny max}}}
%---------------------------------------------------------------------------


%---------------------------------------------------------------------------
\title{Particle filtering for anaerobic bioprocesses}
\author{Fabien Campillo\thanks{MODEMIC INRIA/INRA project-team, Montpellier 
-- \texttt{Fabien.Campillo@inria.fr}}
\and
Jérôme Harmand\thanks{LBE/INRA, Narbonne 
-- \texttt{Jerome.Harmand@supagro.inra.fr}}\ $^*$
\and
Boumédiène Benyahia\thanks{LAT, Université de Tlemcen
-- \texttt{benyahia@supagro.inra.fr}}\ $^*$
}
%---------------------------------------------------------------------------


%---------------------------------------------------------------------------
\maketitle
%---------------------------------------------------------------------------


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{abstract}
%xxxxx
%\paragraph{Keywords and phrases:} 
%xxxx
%\end{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



Particle filtering techniques (sequential Monte Carlo) are called on to develop in the area 
of ​​online monitoring  of biotechnological  process. The robustness of these methods is 
well known and may be useful to deal with noisy measurements, strongly nonlinear dynamics 
but also in case of model mismatch. 

In this preliminary work we present an implementation of the particle filter within the framework of the chemostat.
We perform some tests in the case of simulated observations and when the coefficients are assumed to be known. One of the characteristics of applications in biotechnology is that the observations have a rather low frequency; we focus on that point in the present work.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The problem}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Consider the solution $Z_{t}=(X_{t},S_{t})$ of the two-dimensional stochastic chemostat model \cite{campillo2011chemostat}:
\begin{subequations}
\label{eq.state}
\begin{align}
\label{eq.state.X}
  \rmd X_{t} &= (\mu(S_{t})-D)\,X_{t}\,\rmd t
  +c_1\,\sqrt{X_{t}}\,\rmd W^1_{t}
\\
\label{eq.state.S}
  \rmd S_{t} &= -k\,\mu(S_{t})\,X_{t}\,\rmd t+D\,(\Sin-S_{t})\,\rmd t
  +c_2\,\sqrt{S_{t}}\,\rmd W^2_{t}
\end{align}
\end{subequations}
for $t\in[0,T]$, where $X_{t}$ and $S_{t}$ are the concentrations (g/l) of biomass and substrate at time $t$; $W_{t}^1$ and $W_{t}^2$ are independent scalar standard Brownian motions. The parameters are the dilution rate $D$ (1/h), the input substrate concentration $\Sin$ (mg/l), the stoichiometric coefficient $k$ and the noise intensities $c_{1}$ and $c_{2}$. The initial distribution law of the initial condition $Z_{0}=(X_{0},S_{0})$ is denoted $p_{0}(z)=p_{0}(x,s)=p_{Z_{0}}(z)=p_{X_{0},S_{0}}(x,s)$.


We suppose that the specific growth function is of Monod type:
  $\mu(s) = \mumax\,s/(k_{s}+s) $ with   maximum growth rate $\mumax$ (1/h) and  half-saturation $k_{s}$ (mg/l). 

\bigskip

We have observations at instants $t_{n}=n\,\Delta$
\begin{align}
\label{eq.obs}
  Y_{n} &= S_{t_{n}} + \sigma\,S_{t_{n}} \,v_n
\end{align}
where $v_{n}\simiid N(0,1)$ and $T=\nmax\,\Delta$. The state noises $W^i_{t}$, the observation noise $v_{n}$ and the initial condition $(X_{0},B_{0})$ are supposed independent.


We will use the notations
$f_{1}(z)=f_{1}(x,s)=(\mu(s)-D)\,x$ and
 $f_{2}(x,s)=-k\,\mu(s)\,x+D\,(\Sin-s)$.

\medskip

At time  $t_{n}$, we have measurements  $y_{0},\dots,y_{n}$ (noté $y_{0:n}$) that are realisations of $Y_{0},\dots,Y_{n}$ (denoted $Y_{0:n}$). The filtering problem is to compute the optimal nonlinear optimal filter  $\pi_{n}(z)=\pi_{n}(x,s)$ defined by:
\begin{align}
  \pi_{n}(z) \eqdef p_{Z_{t_{n}}|Y_{0:n}}(z|y_{0:n})
\end{align}
the p.d.f. of $Z_{t_{n}}$ given $Y_{0:n}=y_{0:n}$.
We also defined the predicted filter $
  \pi_{n^-}(z) 
  \eqdef
  p_{Z_{t_{n}}|Y_{0:n-1}}(z|y_{0:n-1})$; $\pi_{n^-}$ (resp. $\pi_{n}$) gathers information available on   $Z_{t_{n}}$ based on the observations $Y_{0:n-1}=y_{0:n-1}$
(resp. $Y_{0:n}=y_{0:n}$).



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The optimal filter and the particle filter}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


The iteration $\pi_{n-1}\to \pi_{n}$ is achieved in two classic steps. First we compute $\pi_{n^-}$ with the \emph{prediction step}:
\begin{align}
\label{eq.nlf.prediction}
  \pi_{n^-}(z')
  &=
  \int_{\R^2_{+}} Q_{\Delta}(z,z')\,\pi_{n-1}(z)\,\rmd z
\end{align}
where  $Q_{\Delta}(z,z')$ is the \emph{transition kernel} of the state equation   \eqref{eq.state} defined by:
\begin{align}
\label{eq.Q}
  Q_{\Delta}(z,z')
  \eqdef
  p_{Z_{t+\Delta}|Z_{t}}(z'|z)
  \,.
\end{align}
Then, using the new observation $Y_{n}=y_{n}$, we compute $\pi_{n}$ with the  \emph{correction step}:
\begin{align}
\label{eq.nlf.correction}
\pi_{n}(z')
  =
  \frac{\Psi(y_{n},z)\,\pi_{n^-}(z)}
        {\int \Psi(y_{n},z')\,\pi_{n^-}(z')\,\rmd z'}
\end{align}
also denoted $  \pi_{n}(z')
  \propto
  \Psi(y_{n},z)\,\pi_{n^-}(z)$, where $\Psi(y,z)$ is the \emph{likelihood function} of the last observation:
\begin{align}
\label{eq.likelihood}
  \Psi(y,z)
  &\eqdef
  \exp\Big(-\frac{1}{2\,\sigma^2\,s^2}\,|y-s|^2\Big)\,.
\end{align}
i.e. the p.d.f. of $Y_{n}=y$ given $Z_{t}=z=(x,s)$. 

The initialization at time $t=0$ of the iterations is $\pi_{0^-}=\pi^0$ and
$  \pi_{0}(z')
  \propto
  \Psi(y_0,z)\,\pi_{0^-}(z)$.


Except in the linear/Gaussian case and some very particular cases, this optimal filter  \eqref{eq.nlf.prediction} and \eqref{eq.nlf.correction} cannot be solve explicitly. This is why specific approximation techniques have been develop. A first method is the extended Kalman filter; another one is the particle filer. We present the latter now.

\bigskip

We introduce the simplest implementation of the particle filter, namely the bootstrap filter \cite{doucet2001a}. The idea is to obtain an empirical approximation of $\pi_{n}$ of the form:
\[
 \pi_{n}(z) 
 \simeq 
 \pi_{n}^N(z) 
 \eqdef 
 \sum_{i=1}^N\delta_{\xi^i_{n}}(z)
\]
where $\xi^1_{n},\dots,\xi^N_{n}$ are the particles, $\delta_{\xi}(z)$ is the Dirac measure centered on  $\xi$. 
The particle filter in a sequential Monte Carlo method where  $\pi_{n}^N$ is computed from  $\pi_{n-1}^N$ and the new observation   $Y_{n}=y_{n}$. 

Starting from the particles $(\xi^i_{n-1})_{i=1\cdots N}$, the new particles   $(\xi^i_{n})_{i=1\cdots N}$ are computed in two steps: the prediction (or mutation step) and the correction (or selection step).

%---------------------------------------------------------------------------
\subsubsection*{\it Prediction (mutation)}
%---------------------------------------------------------------------------


The kernel $Q_{\Delta}$ and the prediction step \eqref{eq.nlf.prediction}
are complex so we replace them by their empirical counterpart.  First we have to simulate predicted particles $\tilde\xi^i_{n} \sim \tilde Q_{\Delta}(\xi^i_{n-1},\cdot)$, but the transition kernel $Q_{\Delta}$ is not explicit and we have to approximate it with an Euler-Maruyama approximation of 
\eqref{eq.state}, that is:
\begin{subequations}
\label{eq.state.euler}
\begin{align}
\label{eq.state.X.euler}
  \tilde X_{\tilde t_{m}} &= 
  \Big[
  \tilde X_{\tilde t_{m-1}}
  +f_{1}(\tilde X_{\tilde t_{m-1}},\tilde S_{\tilde t_{m-1}})\,\deltafilter
  +c_1\,\sqrt{\tilde X_{\tilde t_{m-1}}}\,\sqrt{\deltafilter}\,w^1_{m}
  \Big]_{+}
\\
\label{eq.state.S.euler}
  \tilde S_{\tilde t_{m}} 
  &= 
  \Big[
  \tilde S_{\tilde t_{m-1}} 
  +f_{2}(\tilde X_{\tilde t_{m-1}},\tilde S_{\tilde t_{m-1}})\,\deltafilter
  +c_2\,\sqrt{\tilde S_{\tilde t_{m-1}}}\,\sqrt{\deltafilter}\,w^2_{m}
  \Big]_{+}
\end{align}
\end{subequations}
where $w^1_{n}$ and $w^2_{n}$ are $N(0,1)$ independent random variables;
and $\deltafilter=\Delta/\Mfilter$ and $\tilde t_{m}=m\,\deltafilter$. The positive part operator $[\,\cdot\,]_{+}$ is used to ensure the positivity of the solution. 

Starting from $(\tilde X_{0},\tilde S_{0})=\xi^i_{n-1}$, we simulate $M$ time iterations of \eqref{eq.state.euler} to get the predicted particle $\tilde\xi^i_{n} = (\tilde X_{\tilde t_{M}},\tilde S_{\tilde t_{M}})
= (\tilde X_{\Delta},\tilde S_{\Delta})$. We denote this step by:
\[
   \tilde\xi^i_{n} \sim \tilde Q_{\Delta}(\xi^i_{n-1},\,\cdot\,)\,.
\]
%corresponding to the Algorithm \ref{algo.Q}.


%%--------------------
%\begin{algorithm}
%\begin{center}
%\begin{minipage}{\textwidth}
%\hrulefill\\[-1em]
%%\mbox{}
%\begin{algorithmic}
%\STATE $(x,s)$ given
%\STATE $\deltafilter \ot \Delta/\Mfilter$
%\FOR {$m=1:\Mfilter$}
%  \STATE $w_{1}\sim N(0,1)$, $w_{2}\sim N(0,1)$
%  \STATE $x' \ot \max(0,x+f_{1}(x,s)\,\deltafilter+c_{1}\sqrt{x}\,\sqrt{\deltafilter}\,w_{1})$
%  \STATE $s' \ot \max(0,s+f_{2}(x,s)\,\deltafilter+c_{2}\sqrt{s}\,\sqrt{\deltafilter}\,w_{2})$
%  \STATE $x \ot x'$, $s \ot s'$
%\ENDFOR
%\end{algorithmic}
%\hrulefill
%\end{minipage}
%\end{center}
%\caption{\it Simulation of Equation \eqref{eq.state} with  an Euler-Maruyama scheme leading to the approximation $\tilde Q_{\Delta}$ of the transition kernel $Q_{\Delta}$ given by  \eqref{eq.Q}.}
%\label{algo.Q}
%\end{algorithm}
%%--------------------



%---------------------------------------------------------------------------
\subsection*{\it Correction (selection)}
%---------------------------------------------------------------------------

In a second step a likelihood weight $\omega^i_{n}$ is associated with each predicted particle $\tilde \xi^{i}_{n}$. This weight:
\[
   \omega^i_{n}
   \propto
   \Psi(y_{n},\tilde \xi^{i}_{n})
\]
is proportional to the likelihood of the last observation and $\sum_{i=1}^N\omega^i_{n}=1$. Then these particle are selected accordingly to these weights, i.e. sampled from the empirical distribution 
$\sum_{i=1}^N\omega^i\,\delta_{\tilde\xi^i_{n}}$
independantly one from the other:
\begin{align*}
  \xi^1_{n},\dots,\xi^N_{n}
  \simiid
  \sum_{i=1}^N\omega^i_{n}\,\delta_{\tilde\xi^i_{n}}\,.
\end{align*}




%%--------------------
%\begin{algorithm}
%\begin{center}
%\begin{minipage}{\textwidth}
%\hrulefill\\[-1em]
%%\mbox{}
%\begin{algorithmic}
%\STATE{ }\COMMENT{initialisation}     
%\STATE $\tilde\xi^{1},\dots,\tilde\xi^{N}\simiid \law(Z_{0})$
%\STATE $\omega^i \ot \Psi(y_{0},\tilde\xi^{i})$
%     pour $i=1:N$
%     \COMMENT{poids de vraisemblance}
%     \vphantom{$\int_{0}^N$}
%\STATE $\omega^i \ot \omega^i/\sum_{i'=1}^N\omega^{i'}$
%     pour $i=1:N$
%     \COMMENT{normalisation}
%     \vphantom{$\int_{0}^N$}
%\STATE $\xi^{1},\dots,\xi^{N}\simiid \sum_{i=1}^N\omega^i\,\delta_{\xi^i}$
%     \COMMENT{correction (sélection)}
%     \vphantom{$\int_{0}^N$}
%\STATE sauvegarder $(0,\xi^{1:N})$
%     \vphantom{$\int_{0}^N$}
%\STATE{ }     
%\STATE{ }\COMMENT{itérations}     
%\FOR {$n=1,2,\dots$}
%  \STATE  $\tilde\xi^i \sim \tilde Q_{\Delta}(\xi^i,\rmd z')$
%     pour $i=1:N$
%     \COMMENT{prédiction (mutation)}
%     \vphantom{$\int_{0}^N$}
%  \STATE $\omega^i \ot \Psi(y_{n},\tilde\xi^{i})$
%     pour $i=1:N$
%     \COMMENT{poids de vraisemblance}
%     \vphantom{$\int_{0}^N$}
%  \STATE $\omega^i \ot \omega^i/\sum_{i'=1}^N\omega^{i'}$
%     pour $i=1:N$
%     \COMMENT{normalisation}
%     \vphantom{$\int_{0}^N$}
%  \STATE $\xi^{1},\dots,\xi^{N}\simiid \sum_{i=1}^N\omega^i\,\delta_{\xi^i}$
%     \COMMENT{correction (selection)}
%     \vphantom{$\int_{0}^N$}
%  \STATE sauvegarder $(n\,\Delta,\xi^{1:N})$
%     \vphantom{$\int_{0}^N$}
%\ENDFOR
%\end{algorithmic}
%\hrulefill
%\end{minipage}
%\end{center}
%\caption{\it Filtre particulaire (filtre ``bootstrap''): permet de calculer 
%une approximation de $\pi_{n^-}$ de la forme $\sum_{i=1}^N\delta_{\tilde\xi^i}(\rmd z)$ et une approximation de $\pi_{n}(\rmd z)$ de la forme $\sum_{i=1}^N\delta_{\xi^i}(\rmd z)$.  Dans les étapes de prédiction, les particules prédites $\tilde\xi^i$ sont échantillonnés indépendamment les uns des autres.}
%\label{algo.particle}
%\end{algorithm}
%%--------------------


%\begin{remark}[correction progressive]
%Dans l'Algorithme \ref{algo.particle}, lors\-que la somme des poids $\sum_{i'=1}^N\omega^{i'}$ est nulle cela signifie que le filtre a ``perdu la piste'' de l'état courant. 
%Dans ce cas il existe plusieurs possibilités. En effet, dans la plupart des applications biotechnologiques la fréquence des observations est plutôt faible et il est donc possible de simuler des particules prédites jusqu'à obtenir un ensemble de particules ``suffisamment vraisemblables''. On peut également faire appel à des techniques spécifiques comme celle de la correction progressive, voir 
%\citet[\S\,5.5]{oudjane2000b} et \citet{musso2001a}.
%\end{remark}
%
%
%\begin{remark}[Lissage]
%Si  $T=\nmax\,\Delta$ est le temps final de l'expérimen\-tation ou de l'exploitation du bioréacteur, 
%il peut être intéressant de calculer le lisseur: 
%\[
%  \bar\pi_{n}(z)
%  \eqdef
%  p_{Z_{n}|Y_{0:\nmax}}(z|y_{0:\nmax})\,,
%\]
%voir \cite{doucet2009tutorial,briers2010smoothing}
%\end{remark}






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical test}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In this prelimnray work on simulated observations obtained with an Euler-Maruyama scheme \eqref{eq.state.euler} with a new time step $\deltasimu=\Delta/\Msimu$ where $\Msimu$ is a multiple of $\Mfilter$. Observations are simulated accordingly to \eqref{eq.obs}.






%---------------------------------------------------------------------------
\subsection{Test 1}
%---------------------------------------------------------------------------

Final time is $T=1000$ (h) and we use the same discrtization step for the observations, the simulation of the SDE and the filter: 
$\Delta=\deltasimu=\deltafilter=0.5$ (i.e. $\Msimu=\Mfilter=1$). The parameters are: 
dilution rate $D=0.01\,\textrm{h}^{-1}$; 
input substrate concentration $\Sin=100\,\textrm{mg}/\textrm{l}$; 
stoichiometric coefficient $k=10$; 
maximum growth rate $\mumax=0.3\,\textrm{h}^{-1}$; 
half-saturation  $k_{s}=10\,\textrm{mg}/\textrm{l}$. 
L'instant final  $T=1000$ correspond à 10 fois le temps de rétention $\frac1D$.

The state noise intensities are $c_{1}=c_{2}=0.03$, the observation noise intensity is $\sigma=0.2$. The initial law is $\pi^0(\rmd x,\rmd s)=\NN(0.2,0.5^2)\otimes\NN(1,0.5^2)$.

The number of particles is $N=1000$.






%-------------------
\begin{figure}
\centering
\includegraphics[angle=90,width=6.28cm]{matlab/test1/graph_b_and_s.pdf}
\hfill
\includegraphics[angle=90,width=6cm]{matlab/test1/graph_b_and_s_tubes.pdf}
\caption{Test 1. Two first left: the simulated process $t\to (X_{t},S_{t})$ (blue), the observation process $t_{n}\to Y_{n}$ (green), the estimates $t\to (\hat X_{t},\hat S_{t})$ (red) and the deterministic trajectory $t\to (x(t),s(t))$ (black) 
--- Two last right: same graphics but without the observation process and with the grey ``tubes'' around the estimates $\hat X_{t}$ and $\hat S_{t}$ that are the minimum and maximum values taken by the $x$ and $s$ components of the particles.}
\label{fig.example1a}
\end{figure}
%-------------------




%%-------------------
%\begin{figure}
%\centering
%\includegraphics[width=8cm]{matlab/test1/graph_phase.pdf}
%\caption{Test 1, phase plot output obtained with \texttt{run(4)} (graphic procedure \texttt{graph\symbol{95}phase.pdf}).}
%\label{fig.example1b}
%\end{figure}
%%-------------------


%%-------------------
%\begin{figure}
%\centering
%\includegraphics[width=11cm]{matlab/test1/graph_b_and_s.pdf}
%\caption{Test 1, sorties obtenues avec  \texttt{run(4)}; le processus simulé    $t\to (X_{t},S_{t})$ (bleu), le processus d'observation  $t_{n}\to y_{n}$ (vert), les estimations   $t\to (\hat X_{t},\hat S_{t})$ (rouge) et la trajectoire déterministe  $t\to (x(t),s(t))$ (noir) (procédure graphique \texttt{graph\symbol{95}b\symbol{95}and\symbol{95}s}).}
%\label{fig.example1}
%\end{figure}
%%-------------------
%
%%-------------------
%\begin{figure}
%\centering
%\includegraphics[width=11cm]{matlab/test1/graph_b_and_s_tubes.pdf}
%\caption{Test 1, sorties obtenues à l'aide de \texttt{run(4)}; identique à la  Figure \ref{fig.example1} mais sans le processus d'observation et avec les ``tubes'' de confiance en gris autour des estimations  $\hat X_{t}$ et $\hat S_{t}$ correspondant aux valeurs minimum et maximum prisent par les composantes  $x$ et  $s$ des particules (procédure graphique \texttt{graph\symbol{95}b\symbol{95}and\symbol{95}s\symbol{95}tubes}).}
%\label{fig.example2}
%\end{figure}
%%-------------------
%
%%-------------------
%\begin{figure}
%\centering
%\includegraphics[width=11cm]{matlab/test1/graph_phase.pdf}
%\caption{Test 1, tracé dans l'espace des phases obtenu par  \texttt{run(4)} (procédure graphique \texttt{graph\symbol{95}phase}).}
%\label{fig.example3}
%\end{figure}
%%-------------------


%---------------------------------------------------------------------------
\subsection{Test 2}
%---------------------------------------------------------------------------


This test is identical to Test 1 except that we do not use the same time step for the simulation of the SDE, for the observations and for the filter: $\Delta=10$, $\Delta=\deltasimu=\deltafilter=0.5$ (i.e. $\Msimu=\Mfilter=20$). Note that in this test, the particle filter is much faster than in Test 1, the respective CPU times are 17.3945 (s) for Test 1 and 0.84054 (s) for Test 2. This is du to the fact that the selection/correction step cannot be ``vectorized'' (like the mutation/prediction step) and is very time consuming.





%-------------------
\begin{figure}
\centering
\includegraphics[angle=90,width=5.9cm]{matlab/test2/graph_b_and_s.pdf}
\hfill
\includegraphics[angle=90,width=6.3cm]{matlab/test2/graph_b_and_s_tubes.pdf}
\caption{Test 2. Two first left: 
the simulated process $t\to (X_{t},S_{t})$ (blue), the observation process $t_{n}\to Y_{n}$ (green), the estimates $t\to (\hat X_{t},\hat S_{t})$ (red) and the deterministic trajectory $t\to (x(t),s(t))$ (black) 
--- Two last right: 
same graphics but without the observation process and with the grey ``tubes'' around the estimates $\hat X_{t}$ and $\hat S_{t}$ that are the minimum and maximum values taken by the $x$ and $s$ components of the particles.}
\label{fig.example2a}
\end{figure}
%-------------------



%%-------------------
%\begin{figure}
%\centering
%\includegraphics[width=8cm]{matlab/test2/graph_phase.pdf}
%\caption{Test 2, phase plot output obtained with \texttt{run(4)}.}
%\label{fig.test2.example2a}
%\end{figure}
%%-------------------




%%-------------------
%\begin{figure}
%\centering
%\includegraphics[width=11cm]{matlab/test2/graph_b_and_s.pdf}
%\caption{Test 2, sorties obtenues avec \texttt{run(4)}; processus simulé $t\to (X_{t},S_{t})$ (bleu), processus d'observation $t_{n}\to y_{n}$ (vert), les estimations  $t\to (\hat X_{t},\hat S_{t})$ (rouge) et la trajectoire déterministe  $t\to (x(t),s(t))$ (noir) (procédure graphique \texttt{graph\symbol{95}b\symbol{95}and\symbol{95}s}).}
%\label{fig.test2.example1}
%\end{figure}
%%-------------------
%
%%-------------------
%\begin{figure}
%\centering
%\includegraphics[width=11cm]{matlab/test2/graph_b_and_s_tubes.pdf}
%\caption{Test 2, sorties obtenues avec \texttt{run(4)}; identique à la  Figure \ref{fig.test2.example1} mais sans le processus d'observation et avec les ``tubes'' de confiance en gris autour des estimations  $\hat X_{t}$ et $\hat S_{t}$ correspondant aux valeurs minimum et maximum prisent par les composantes  $x$ et  $s$ des particules (procédure graphique \texttt{graph\symbol{95}b\symbol{95}and\symbol{95}s}).}
%\label{fig.test2.example2}
%\end{figure}
%%-------------------
%
%%-------------------
%\begin{figure}
%\centering
%\includegraphics[width=11cm]{matlab/test2/graph_phase.pdf}
%\caption{Test 2, tracé dans l'espace des phases obtenu par  \texttt{run(4)} (procédure graphique \texttt{graph\symbol{95}phase}).}
%\label{fig.test2.example3}
%\end{figure}
%%-------------------
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Discussion and perspectives}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We have developed a first implementation of a particle filter in a chemostat model. One challenge is to manage the fact that the frequency of observations is low, hence it is necessary to simulate the state equation  between two observations. This prediction step can be done in an efficient way as it is highly parallelizable.


Perspectives are numerous. The first is that it is necessary to develop a joint procedure for filtering and parameter estimation. The second is that it would be pertinent to calculate the smoother, i.e. the p.d.f. of $Z_{t_{n}}$ given all the observations $y_{0},\dots,y_{\nmax}$; the corresponding algorithm is more complex. The most interesting perspective is to test the robustness of the filter in case of miss match model: use a filter with a given growth function $\mu(s)$ as the observations correspond to a different growth function $\tilde\mu(s)$. Another interesting question is  to determine which growth function, among a given finite set of growth functions $\{\mu_{i}(s);i=1,\dots,I\}$, is underlying a given set of of observations.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\nocite{venkateswarlu2004a}
\addcontentsline{toc}{section}{Reference}
\bibliographystyle{plainnat}
\bibliography{lib/fab}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



